*
* $Id$
*
* $Log: gdrelx.F,v $
* Revision 1.1.1.1  2002/07/24 15:56:25  rdm
* initial import into CVS
*
* Revision 1.1.1.1  2002/06/16 15:18:41  hristov
* Separate distribution  of Geant3
*
* Revision 1.1.1.1  1999/05/18 15:55:20  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:21:24  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/03 06/10/94  18.22.43  by  S.Giani
*-- Author :
      SUBROUTINE G3DRELX(A,Z,DENS,T,HMASS,DEDX)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  Calculates the mean 1/DENS*dE/dx of a particle with kinetic   *
C.    *  energy T in an element of atomic number Z, atomic weight A    *
C.    *  and density DENS ( the density is just used for the           *
C.    *  calculation of the density effect in the case of high T).     *
C.    *  The routine reproduces the experimental and/or tabulated      *
C.    *  energy losses rather well down to T -> 0.                     *
C.    *  Simple parametrization is used for  T .le. T2L=2 MeV (see     *
C.    *  H.H.Andersen,J.F.Ziegler:Hydrogen stopping powers and         *
C.    *  ranges in all element,Pergamon Press,1977.).                  *
C.    *  For T .gt. T2L=2 MeV the corrected Bethe-Bloch stopping       *
C.    *  power / restricted energy loss formula is used.               *
C.    *                                                                *
C.    *                                                                *
C.    *    ==>Called by : G3DRELA                                      *
C.    *       Author    L.Urban    *********                           *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gconsp.inc"
#include "geant321/gccuts.inc"
#include "geant321/gcunit.inc"
      PARAMETER (AMUKEV=931494.32,D=0.000153537,T1L=0.00001,T2L=0.002)
      PARAMETER (EMPROT=0.9382723)
      DIMENSION B(6,92),C(6,92),CECOF(6)
*
      DATA ((B(I,J),I=1,6),J=1,19) /
     +     1.262,1.44,242.6,12000.,0.1159,18.8,
     +     1.229,1.397,484.5,5873.,0.05225,41.7,
     +     1.411,1.6,725.6,3013.,0.04578,47.6,
     +     2.248,2.59,966.,153.8,0.03475,62.7,
     +     2.474,2.815,1206.,1060.,0.02855,76.0,
     +     2.631,2.989,1445.,957.2,0.02819,77.3,
     +     2.954,3.35,1683.,1900.,0.02513,86.7,
     +     2.652,3.,1920.,2000.,0.0223,97.7,
     +     2.085,2.352,2157.,2634.,0.01816,120.,
     +     1.951,2.199,2393.,2699.,0.01568,139.,
     +     2.542,2.869,2628.,1854.,0.01472,148.,
     +     3.792,4.293,2862.,1009.,0.01397,156.,
     +     4.154,4.739,2766.,164.5,0.02023,162.,
     +     4.15,4.7,3329.,550.,0.01321,165.,
     +     3.232,3.647,3561.,1560.,0.01267,172.,
     +     3.447,3.891,3792.,1219.,0.01211,180.,
     +     5.047,5.714,4023.,878.6,0.01178,185.,
     +     5.731,6.5,4253.,530.,0.01123,194.,
     +     5.151,5.833,4482.,545.7,0.01129,193./
      DATA ((B(I,J),I=1,6),J=20,38) /
     +     5.521,6.252,4710.,553.3,0.01112,196.,
     +     5.201,5.884,4938.,560.9,0.009995,218.,
     +     4.862,5.496,5165.,568.5,0.009474,230.,
     +     4.48,5.055,5391.,952.3,0.009117,239.,
     +     3.983,4.489,5616.,1336.,0.008413,259.,
     +     3.469,3.907,5725.,1461.,0.008829,270.,
     +     3.519,3.963,6065.,1243.,0.007782,280.,
     +     3.14,3.535,6288.,1372.,0.007361,296.,
     +     3.553,4.004,6205.,555.1,0.008763,310.,
     +     3.696,4.175,4673.,387.8,0.02188,322.,
     +     4.21,4.75,6953.,295.2,0.006809,320.,
     +     5.041,5.697,7173.,202.6,0.006725,324.,
     +     5.554,6.3,6496.,110.,0.009689,330.,
     +     5.323,6.012,7611.,292.5,0.006447,338.,
     +     5.874,6.656,7395.,117.5,0.007684,340.,
     +     5.611,6.335,8046.,365.2,0.006244,349.,
     +     6.411,7.25,8262.,220.,0.006087,358.,
     +     5.694,6.429,8478.,292.9,0.006087,358.,
     +     6.339,7.159,8693.,330.3,0.006003,363./
      DATA ((B(I,J),I=1,6),J=39,57)  /
     +     6.407,7.234,8907.,367.8,0.005889,370.,
     +     6.734,7.603,9120.,405.2,0.005765,378.,
     +     6.902,7.791,9333.,442.7,0.005587,390.,
     +     6.425,7.248,9545.,480.2,0.005367,406.,
     +     6.799,7.671,9756.,517.6,0.005315,410.,
     +     6.108,6.887,9966.,555.1,0.005151,423.,
     +     5.924,6.677,10180.,592.5,0.004919,443.,
     +     5.238,5.9,10380.,630.,0.004758,458.,
     +     5.623,6.354,7160.,337.6,0.01394,466.,
     +     5.814,6.554,10800.,355.5,0.004626,471.,
     +     6.23,7.024,11010.,370.9,0.00454,480.,
     +     6.41,7.227,11210.,386.4,0.004474,487.,
     +     7.5,8.48,8608.,348.,0.009074,494.,
     +     6.979,7.871,11620.,392.4,0.004402,495.,
     +     7.725,8.716,11830.,394.8,0.004376,498.,
     +     8.231,9.289,12030.,397.3,0.004384,497.,
     +     7.287,8.218,12230.,399.7,0.004447,490.,
     +     7.899,8.911,12430.,402.1,0.004511,483.,
     +     8.041,9.071,12630.,404.5,0.00454,480./
      DATA ((B(I,J),I=1,6),J=58,76)  /
     +     7.489,8.444,12830.,406.9,0.00442,493.,
     +     7.291,8.219,13030.,409.3,0.004298,507.,
     +     7.098,8.,13230.,411.8,0.004182,521.,
     +     6.91,7.786,13430.,414.2,0.004058,537.,
     +     6.728,7.58,13620.,416.6,0.003976,548.,
     +     6.551,7.38,13820.,419.,0.003877,562.,
     +     6.739,7.592,14020.,421.4,0.003863,564.,
     +     6.212,6.996,14210.,423.9,0.003725,585.,
     +     5.517,6.21,14400.,426.3,0.003632,600.,
     +     5.219,5.874,14600.,428.7,0.003498,623.,
     +     5.071,5.706,14790.,433.,0.003405,640.,
     +     4.926,5.542,14980.,433.5,0.003342,652.,
     +     4.787,5.386,15170.,435.9,0.003292,662.,
     +     4.893,5.505,15360.,438.4,0.003243,672.,
     +     5.028,5.657,15550.,440.8,0.003195,682.,
     +     4.738,5.329,15740.,443.2,0.003186,684.,
     +     4.574,5.144,15930.,442.4,0.003144,693.,
     +     5.2,5.851,16120.,441.6,0.003122,698.,
     +     5.07,5.704,16300.,440.9,0.003082,707./
      DATA ((B(I,J),I=1,6),J=77,92)  /
     +     4.945,5.563,16490.,440.1,0.002965,735.,
     +     4.476,5.034,16670.,439.3,0.002871,759.,
     +     4.856,5.46,18320.,438.5,0.002542,755.,
     +     4.308,4.843,17040.,487.8,0.002882,756.,
     +     4.723,5.311,17220.,537.,0.002913,748.,
     +     5.319,5.982,17400.,586.3,0.002871,759.,
     +     5.956,6.7,17800.,677.,0.00266,765.,
     +     6.158,6.928,17770.,586.3,0.002812,775.,
     +     6.204,6.979,17950.,586.3,0.002776,785.,
     +     6.181,6.954,18120.,586.3,0.002748,793.,
     +     6.949,7.82,18300.,586.3,0.002737,796.,
     +     7.506,8.448,18480.,586.3,0.002727,799.,
     +     7.649,8.609,18660.,586.3,0.002697,808.,
     +     7.71,8.679,18830.,586.3,0.002641,825.,
     +     7.407,8.336,19010.,586.3,0.002603,837.,
     +     7.29,8.204,19180.,586.3,0.002573,847./
      DATA C/92*0.,92*0.,92*0.,92*0.,92*0.,92*0./
      DATA CECOF/0.42237,0.0304,-0.00038,3.858,-0.1668,0.00158/
      DATA EPS/0.000001/
*     ------------------------------------------------------------------
*      in the case of non-integer Z the low energy parameters
*      and the ionization potential are taken at INT(Z) !
*
      IZ=INT(Z+EPS)
      IF((IZ.LT.1).OR.(IZ.GT.92)) GOTO 10
*
*     Calculate coefficients C(I,J) if it has not been done already
*
      IF(C(1,IZ).LT.EPS) THEN
         FAC=AVO/A
         C(1,IZ)=FAC*AMUKEV**0.5*B(1,IZ)
         C(2,IZ)=FAC*AMUKEV**0.45*B(2,IZ)
         C(3,IZ)=FAC*B(3,IZ)/AMUKEV
         C(4,IZ)=B(4,IZ)/AMUKEV
         C(5,IZ)=AMUKEV*B(5,IZ)
*
*                     POTI=16.E-9*Z**0.9
*
         POTI=B(6,IZ)*1.E-9
*
*
         C(6,IZ)=POTI
*
      ENDIF
*     ----------------------------------------------------------------
      T1LIM=HMASS*T1L/EMPROT
      T2LIM=HMASS*T2L/EMPROT
*
*     Calculate dE/dx
*     ---> for T .le. T1LIM (very low energy)
*
      IF(T.LE.T1LIM) THEN
         TAU=T/HMASS
         DEDX=C(1,IZ)*TAU**0.5
      ELSE
*
*     ---> for T1LIM .lt. T   and  T .le. T2LIM (low energy)
*
         IF(T.LE.T2LIM) THEN
            TAU=T/HMASS
            SL=C(2,IZ)*TAU**0.45
            SH=C(3,IZ)*LOG(1.+C(4,IZ)/TAU+C(5,IZ)*TAU)/TAU
            DEDX=SL*SH/(SL+SH)
*
*     ---> for T .gt. T2LIM ( "high " energy , Bethe-Bloch formula)
*
         ELSE
            P=SQRT(T*(T+2.*HMASS))
            E=T+HMASS
            BET2=(P/E)**2
            ETA=P/HMASS
            ETA2=ETA*ETA
*+++ new line follows.....
            B2G2=ETA*ETA
*+++ end of correction
            TMAX=2.*EMASS*T*(T+2.*HMASS)
*+++  correction of the next line
*           TMAX=TMAX/(HMASS**2+EMASS**2+EMASS*(T+HMASS))
            TMAX=TMAX/(HMASS**2+EMASS**2+2.*EMASS*E)
*+++ end of correction
*
*         density correction
*
            POTI=C(6,IZ)
            CC=1.+2.*LOG(POTI/(28.8E-9*SQRT(DENS*Z/A)))
*         condensed material ? ( dens .gt. 0.05 ? )
            IF(DENS.GT.0.05) THEN
               IF(POTI.LT.1.E-7) THEN
                  IF(CC.LT.3.681) THEN
                     X0=0.2
                  ELSE
                     X0=0.326*CC-1.
                  ENDIF
                  X1=2.
               ELSE
                  IF(CC.LT.5.215) THEN
                     X0=0.2
                  ELSE
                     X0=0.326*CC-1.5
                  ENDIF
                  X1=3.
               ENDIF
*         gas ?   ( dens . le . 0.05 ? )
            ELSE
               IF(CC.LE.12.25) THEN
                  IP=INT((CC-10.)/0.5)+1
                  IF(IP.LT.0) IP=0
                  IF(IP.GT.4) IP=4
                  X0=1.6+0.1*FLOAT(IP)
                  X1=4.
               ELSE
                  IF(CC.LE.13.804) THEN
                     X0=2.
                     X1=5.
                  ELSE
                     X0=0.326*CC-2.5
                     X1=5.
                  ENDIF
               ENDIF
            ENDIF
*
            XA=CC/4.606
            XM=3.
            AA=4.606*(XA-X0)/(X1-X0)**XM
*
            X=LOG10(ETA)
            DELTA=0.
            IF(X.GT.X0) THEN
               DELTA=4.606*X-CC
               IF(X.LT.X1) DELTA=DELTA+AA*(X1-X)**XM
            ENDIF
*
*         calculate shell correction
*
            POTSQ=POTI*POTI
            IF(ETA.GT.0.13) THEN
               F1=1./ETA2
               F2=F1*F1
               F3=F1*F2
               F4=(F1*CECOF(1)+F2*CECOF(2)+F3*CECOF(3))*1.E+12
               F5=(F1*CECOF(4)+F2*CECOF(5)+F3*CECOF(6))*1.E+18
               CE=F4*POTSQ+F5*POTSQ*POTI
            ELSE
               ETA2=0.0169
               F1=1./ETA2
               F2=F1*F1
               F3=F1*F2
               F4=(F1*CECOF(1)+F2*CECOF(2)+F3*CECOF(3))*1.E+12
               F5=(F1*CECOF(4)+F2*CECOF(5)+F3*CECOF(6))*1.E+18
               CE=F4*POTSQ+F5*POTSQ*POTI
               CE=CE*LOG(T/T2LIM)/LOG(0.0079/T2LIM)
            ENDIF
*
            F1=D*Z/(A*BET2)
*
*         stopping power or restricted dE/dx ?
*
*+++  correction of the next few lines
*           IF(DCUTM.GE.TMAX) THEN
*              F2=2.*(LOG(TMAX/POTI)-BET2)
*           ELSE
*              F2=LOG(TMAX*DCUTM/POTSQ)-BET2*(1.+DCUTM/TMAX)
*           ENDIF
            TUPP=DCUTM
            IF(TMAX.LT.DCUTM) TUPP=TMAX
            F2=LOG(2.*EMASS*B2G2/POTI)+LOG(TUPP/POTI)
     +         -BET2*(1.+TUPP/TMAX)
*+++ end of correction
            DEDX=F1*(F2-DELTA-2.*CE/Z)
*
*
            TAU=T2LIM/HMASS
            SL=C(2,IZ)*TAU**0.45
            SH=C(3,IZ)*LOG(1.+C(4,IZ)/TAU+C(5,IZ)*TAU)/TAU
            ST=SL*SH/(SL+SH)
*
            TMAX=2.*EMASS*T2LIM*(T2LIM+2.*HMASS)
*+++  correction of the next line
*           TMAX=TMAX/(HMASS**2+EMASS**2+EMASS*(T2LIM+HMASS))
            TMAX=TMAX/(HMASS**2+EMASS**2+2.*EMASS*(T2LIM+HMASS))
*+++  end of correction
            BET2=T2LIM*(T2LIM+2.*HMASS)/(T2LIM+HMASS)**2
            SBB=2.*(LOG(TMAX/POTI)-BET2)
            SBB=D*Z*SBB/(A*BET2)
            CORBB=(ST/SBB-1.)*T2LIM
*
            DEDX=DEDX*(1.+CORBB/T)
*
         ENDIF
      ENDIF
*
      RETURN
   10 WRITE (CHMAIL,10000) Z
      CALL GMAIL(1,1)
10000 FORMAT(5X,'***GDRELP  (Z.LT.1.).OR.(Z.GT.92.) ==> Z=',F10.1)
      END
